/*

  Simple cross-product, outputs in correct triangular form.

  - Coalescing memory access in the Column reads
  - No memory coalescing in the Row reads
  - 50% coalescing in the write (float8s)
  - All threads read in all data

 */

CUBE_KERNEL(shared1x1float4, cuFloatComplex *array_d, cuFloatComplex *product_d,
	    const int write)
{
  CUBE_START;

  //set frequency number from blockIdx.y
  int f = blockIdx.y;

  int Row, Col, blockX, blockY;
  CUBE_DEVICE_CALL(findPosition, Col, Row, blockX, blockY);

  //instantiate sum variables
  float sumXXreal = 0.0, sumXXimag = 0.0;
  float sumXYreal = 0.0, sumXYimag = 0.0;
  float sumYXreal = 0.0, sumYXimag = 0.0;
  float sumYYreal = 0.0, sumYYimag = 0.0;

  //get local thread ID
  int ty = threadIdx.y;
  int tx = threadIdx.x;

  //declare shared memory for input coalescing
  __shared__ float input[TILE_WIDTH*8];

  //  #pragma unroll
  for(int t=0; t<NTIME; t++){
    __syncthreads();
    int tid = ty*TILE_WIDTH + tx;
    if(tid < TILE_WIDTH){
      float4 tempCol = ((float4*)array_d)[t*NSTATION*NFREQUENCY + f*NSTATION + Col];
      input[tx] = tempCol.x;
      input[tx + TILE_WIDTH] = tempCol.y;
      input[tx + 2*TILE_WIDTH] = tempCol.z;
      input[tx + 3*TILE_WIDTH] = tempCol.w;
      CUBE_ADD_BYTES(16);
    }
    if(tid >= TILE_WIDTH && tid < (TILE_HEIGHT + TILE_WIDTH)){
      float4 tempRow = ((float4*)array_d)[t*NSTATION*NFREQUENCY + f*NSTATION + blockY*TILE_HEIGHT + tx];
      input[tx + TILE_WIDTH*4] = tempRow.x;
      input[tx + TILE_HEIGHT + TILE_WIDTH*4] = tempRow.y;
      input[tx + 2*TILE_HEIGHT + TILE_WIDTH*4] = tempRow.z;
      input[tx + 3*TILE_HEIGHT + TILE_WIDTH*4] = tempRow.w;
      CUBE_ADD_BYTES(16);
    }
    __syncthreads();

    float colXreal = input[tx];
    float colXimag = input[tx + 1*TILE_WIDTH];
    float colYreal = input[tx + 2*TILE_WIDTH];
    float colYimag = input[tx + 3*TILE_WIDTH];
    float rowXreal = input[ty + 4*TILE_WIDTH];
    float rowXimag = input[ty + 5*TILE_WIDTH];
    float rowYreal = input[ty + 6*TILE_WIDTH];
    float rowYimag = input[ty + 7*TILE_WIDTH];

    sumXXreal += rowXreal * colXreal;
    sumXXreal += rowXimag * colXimag;
    sumXXimag += rowXimag * colXreal;
    sumXXimag -= rowXreal * colXimag;
    sumXYreal += rowXreal * colYreal;
    sumXYreal += rowXimag * colYimag;
    sumXYimag += rowXimag * colYreal;
    sumXYimag -= rowXreal * colYimag;
    sumYXreal += rowYreal * colXreal;
    sumYXreal += rowYimag * colXimag;
    sumYXimag += rowYimag * colXreal;
    sumYXimag -= rowYreal * colXimag;
    sumYYreal += rowYreal * colYreal;
    sumYYreal += rowYimag * colYimag;
    sumYYimag += rowYimag * colYreal;
    sumYYimag -= rowYreal * colYimag;

    /*    sumXXreal += rowXreal * colXreal + rowXimag * colXimag;
    sumXXimag += rowXimag * colXreal - rowXreal * colXimag;
    sumXYreal += rowXreal * colYreal + rowXimag * colYimag;
    sumXYimag += rowXimag * colYreal - rowXreal * colYimag;
    sumYXreal += rowYreal * colXreal + rowYimag * colXimag;
    sumYXimag += rowYimag * colXreal - rowYreal * colXimag;
    sumYYreal += rowYreal * colYreal + rowYimag * colYimag;
    sumYYimag += rowYimag * colYreal - rowYreal * colYimag;*/
  } 

  if (write)
    CUBE_DEVICE_CALL(write1x1, Col, Row, product_d,  
		     sumXXreal,  sumXXimag, sumXYreal,  sumXYimag,  
		     sumYXreal, sumYXimag,  sumYYreal,  sumYYimag);

  if (Col <= Row) {
    CUBE_ADD_BYTES(32);
    CUBE_ADD_FLOPS(NTIME*32);
  }

  CUBE_END;
}
